 #This script performs OLS and quantile (linear) regressions for 4 species of skinks, and prepares a datafile to compare residual abundance to microclimate

#Clear workspace

rm(list=ls())

in.dir = "C:/R/In/"
out.dir = "C:/R/Out/"

#Set working directory
setwd(in.dir)

#load necessary libraries

#list the libraries needed
necessary=c("adehabitat","SDMTools","quantreg")
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)

#load the libraries
for (lib in necessary) library(lib,character.only=T)

#Import the dataset for analysis

rawdata= read.csv("pairedabundanceandenvironmentalsuitabilitybysite.csv",header=T)

#Use thresholds reported in MAXENT to remove data where the ES value is below the threshold, start by importing threshold file

max.thresh = read.csv("Species Threshold Values from MAXENT.csv", header=T)

#rawdata$carrubr_es[which(rawdata$carrubr_es<max.thresh[3,2])] = NA 
#rawdata$lamcogg_es[which(rawdata$lamcogg_es<max.thresh[1,2])] = NA
#rawdata$sapbasi_es[which(rawdata$sapbasi_es<max.thresh[2,2])] = NA
#rawdata$gnyquee_es[which(rawdata$gnyquee_es<max.thresh[4,2])] = NA

#Perform OLS (Ordinary Least Squares Regression) on dataset, comparing each species ES value to abundance and SD of abundance

lcESabund = lm(lamcogg_abund~lamcogg_es, data=rawdata)
gqESabund = lm(gnyquee_abund~gnyquee_es, data=rawdata)
crESabund = lm(carrubr_abund~carrubr_es, data=rawdata)
sbESabund = lm(sapbasi_abund~sapbasi_es, data=rawdata)

lcESsd = lm(lamcogg_sd~lamcogg_es, data=rawdata)
gqESsd = lm(gnyquee_sd~gnyquee_es, data=rawdata)
crESsd = lm(carrubr_sd~carrubr_es, data=rawdata)
sbESsd = lm(sapbasi_sd~sapbasi_es, data=rawdata)

#Prepare a summary table of OLS coefficients

OLSsummary = data.frame(spp=c("LAMCOGG","GNYQUEE","CARRUBR","SAPBASI"), int=c(coefficients(lcESabund)[1],coefficients(gqESabund)[1], coefficients(crESabund)[1], coefficients(sbESabund)[1]), slope=c(coefficients(lcESabund)[2],coefficients(gqESabund)[2], coefficients(crESabund)[2], coefficients(sbESabund)[2]), Rsqr=c(summary(lcESabund)$r.squared, summary(gqESabund)$r.squared, summary(crESabund)$r.squared, summary(sbESabund)$r.squared))

#Perform quantile regression on dataset, comparing each species ES value to abundance and SD

QlcESabund = rq(lamcogg_abund~lamcogg_es, tau=c(.95), data=rawdata, method="br", model=T)
QgqESabund = rq(gnyquee_abund~gnyquee_es, tau=c(.95), data=rawdata, method="br", model=T)
QcrESabund = rq(carrubr_abund~carrubr_es, tau=c(.95), data=rawdata, method="br", model=T)
QsbESabund = rq(sapbasi_abund~sapbasi_es, tau=c(.95), data=rawdata, method="br", model=T)

QlcESsd = rq(lamcogg_sd~lamcogg_es, tau=c(.95), data=rawdata, method="br", model=T)
QgqESsd = rq(gnyquee_sd~gnyquee_es, tau=c(.95), data=rawdata, method="br", model=T)
QcrESsd = rq(carrubr_sd~carrubr_es, tau=c(.95), data=rawdata, method="br", model=T)
QsbESsd = rq(sapbasi_sd~sapbasi_es, tau=c(.95), data=rawdata, method="br", model=T)

#Produce plots of relationships to check with R-viewer

png("C:/R/Out/delete/test2.png", height=800, width=800)
par(mfrow=c(2,2))
plot(rawdata$lamcogg_es, rawdata$lamcogg_abund, xlab="es", ylab="abund", main="LAMCOGG",abline(lcESabund, col="red"))
abline(QlcESabund, col="blue")
plot(rawdata$gnyquee_es, rawdata$gnyquee_abund, xlab="es", ylab="abund", main="GNYQUEE",abline(gqESabund, col="red"))
abline(QgqESabund, col="blue")
plot(rawdata$carrubr_es, rawdata$carrubr_abund, xlab="es", ylab="abund", main="CARRUBR",abline(crESabund, col="red"))
abline(QcrESabund, col="blue")
plot(rawdata$sapbasi_es, rawdata$sapbasi_abund, xlab="es", ylab="abund", main="SAPBASI",abline(sbESabund, col="red"))
abline(QsbESabund, col="blue")
dev.off()

#Create summary table for quantile regressions ***First learn to calculate Neagles pseudo-rsquared

#Create a dataframe with georef_ID, east, north, measured abundance, residual abundance.  Write this dataframe out as a .csv file to examine it's relationships

abund.data = data.frame(georef_ID=rawdata$georef_ID,east=rawdata$east, north=rawdata$north, lcSD=rawdata$lamcogg_sd, lclocalabund=rawdata$lamcogg_abund, gqSD=rawdata$gnyquee_sd, gqlocalabund=rawdata$gnyquee_abund, crSD=rawdata$carrubr_sd, crlocalabund=rawdata$carrubr_abund, sbSD=rawdata$sapbasi_sd, sblocalabund=rawdata$sapbasi_abund)
lcresid.data = data.frame(georef_ID=rawdata$georef_ID, lcresidabund=QlcESabund$residuals)
gqresid.data = data.frame(georef_ID=rawdata$georef_ID, gqresidabund=QgqESabund$residuals)
crresid.data = data.frame(georef_ID=rawdata$georef_ID, crresidabund=QcrESabund$residuals)
sbresid.data = data.frame(georef_ID=rawdata$georef_ID, sbresidabund=QsbESabund$residuals)
data.out=merge(abund.data,lcresid.data,by=c("georef_ID"))
data.out2=merge(data.out,gqresid.data, by=c("georef_ID"))
data.out3=merge(data.out2, crresid.data, by=c("georef_ID"))
data.out4=merge(data.out3,sbresid.data, by=c("georef_ID"))


write.csv(x=data.out4,file=paste(out.dir,"residualabundanceplussd.csv", sep=""), row.names=F)
